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1.0  Introduction 

The  problem  which  we  will  be  considering  is  that  of  interpolation  in 
two  (or  more)  variables,  where  the  data  to  be  interpolated  is  at  irregularly 
spaced  points.   While  the  discussion  will  often  center  on  functions  of  two 
variables,  most  of  the  ideas  are  readily  extendable  to  any  number  of  inde- 
pendent  variables. 

Assume  that  data  points   (X,  ,  f,  )  ,  k  =  1,...,N  are  specified,  where 
the  X,   lie  in  n-space.   We  make  no  assumption  about  the  disposition  of 
the  points  X,  .  We  wish  to  construct  a  function  F(X)  which  is  defined 
in  some  region  which  contains  all  the  X   and  such  that  F(X,  )  ■  f,  , 
k  =  1,...,N  .   We  wish  to  discuss  a  broad  enough  class  of  interpolating 
functions  so  that  we  can  require  that  F(X)   is  not  only  continuous,  but 
perhaps  has  some  continuous  derivatives  as  well. 

A  number  of  methods  have  been  proposed  for  the  problem.   For  the  most 
part  these  are  global  methods,  most  of  which  require  the  solution  of  a  sys- 
tem of  N  linear  equations  for  the  coefficients  of  a  set  of  basis  functions. 
Typical  of  these  are  polynomial  interpolation  [12],  and  optimal  approximations 
in  various  spaces  of  functions   [6],  [7],  [9],  [2], 

Another  type  of  global  approximation  (although  the  computational  version 
is  local)  is  the  idea  of  taking  the  function  value  to  be  a  weighted  average 

of  the  values  at  data  points,  the  weight  being  a  decreasing  function  of  dis- 

-2 
tance  [11].   A  typical  weight  function  is  w(d)  =  d  '  .   There  are  several 

serious  faults  with  this  approximation,  such  as  the  first  partial  derivative 

being  zero  at  each  data  point.   These  faults  are  overcome  by  the  imposition 

of  somewhat  artificial  conditions  on  a  variation  of  the  approximation.   A 


remaining  bad  feature  is  that  the  approximation  is  dependent  on  the  coordi- 
nate system,  or  rather,  the  measure  of  distance. 

A  variation  of  the  distance  weighting  idea  is  proposed  by  McLain  [5]. 
There  the  idea  is  to  fit  the  data  with  an  equation  of  specified  type  in  the 
weighted  least  squares  sense.   The  weight  attached  to  each  point  is  a  de- 
creasing function  of  its  distance  from  the  point  at  which  the  approximation 
is  desired.   McLain  investigated  a  number  of  different  fitting  functions  and 
weight  functions.   This  approach  eliminates  in  a  natural  way  most  of  the  bad 
features  of  Shepard's  approximation,  although  not  the  dependence  on  the  co- 
ordinate system.   McLain  claims  the  approximation  has  partial  derivatives  of 
all  orders,  and  although  this  seems  plausible,  it  is  neither  obvious  nor  eas- 
ily verifiable. 

Where  the  data  is  on  a  somewhat  regular  set  of  points,  e.g.,  if  the  grid 
is  topologically  equivalent  to  a  regular  grid  [3],  or  if  the  grid  can  be  tri- 
angulated, or  represented  as  a  union  of  convex  quadrilaterals  [4],  special 
approaches  may  give  good  results.   Approximations  in  these  cases  would  gen- 
erally be  local  with  limited  smoothness,  although  Ferguson's  [3]  parametric 
representation  is  global  and  a  type  of  spline. 

Another  approach  to  the  interpolation  problem,  the  idea  of  which  is  the 
basis  for  most  of  our  investigation,  is  due  to  Maude  [8].   The  resulting  ap- 
proximation can  be  made  to  have  arbitrarily  high  smoothness  and  is  locally 
determined.   Since  we  will  discuss  the  idea  at  length  in  the  next  section, 
it  is  mentioned  here  only  for  completeness. 

The  ideas  which  will  be  important  to  our  discussion  will  include  the 
following:   (i)  No  assumption  will  be  made  about  the  disposition  of  the  data 
points,  and  in  particular  we  do  not  preclude  the  possibility  of  a  regular 


grid  of  points;  (ii)  A  smooth  approximation  in  the  sense  that  the  function 
and  some  of  its  partial  derivatives  are  continuous ,  is  desirable;  (iii)  The 
number  of  points,   N  ,  may  be  so  large  that  the  approximation  must  be  lo- 
cally determined,  or  at  least  have  a  local  basis  (such  as  the  B-splines  are 
a  local  basis  for  univariate  spline  approximation)  to  be  computationally 
viable;  (iv)  It  is  desirable  for  the  approximation  to  be  independent  of 
translation  and  contraction  or  expansion  of  the  coordinate  system;  and  (v) 
It  is  desirable  that  the  procedure  be  robust  so  that  failure  of  the  exis- 
tence of  the  approximation  can  perhaps  be  circumvented  by  some  slight  al- 
teration. 

It  is  likely  that  (iii)  rules  out  the  optimal  approximations  because 
they  are  global,  and  the  basis  functions  are  sufficiently  complex  that  it 
is  quite  unlikely  a  set  of  basis  functions  can  be  derived  which  have  com- 
pact support.   A  large  number  of  points  also  rules  out  polynomial  approxi- 
mation. 

2.0  A  Class  of  Interpolation  Methods 

We  will  begin  this  section  with  a  brief  review  of  a  method  due  to  Maude 
[8].   Some  of  the  difficulties  experienced  with  the  method  will  be  mentioned. 
A  general  class  of  methods  which  contains  Maude's  will  be  developed  and  some 
properties  of  the  class  of  methods  will  be  discussed. 

2.1  Maude's  Method 

Let   (x,  ,  y.)   and  f  =  f(x,  ,  y,)  ,  k  =  1,...,N  be  specified.   With 

each  point   (x  ,  y  )   associate  a  circle  C   ,  of  radius  r   with  center  at 
n   n  n  n 

(x  ,  y  )  .   This  circle  is  chosen  so  that  the  closed  disk  contains   (x  ,y  ) 
n  Jn  n  n 

and  its  five  nearest  neighbors  from  the  set  {(x,  ,  y,  )}  ,  with  no  more  than 
four  of  the  neighbors  in  the  interior.   There  may  be  'ties'  for  the  position 


of  fifth  nearest  point,  and  one  needs  some  means  of  deciding  which  point  on 
the  boundary  is  to  be  used  in  the  subsequent  calculations.  Let  the  second 
degree  polynomial  interpolating  the  function  f(x,y)   at   (x  ,y  )   and  its 
five  nearest  neighbors  be  denoted  by  Q  (x,y)  .   Let 

12     3 
l-3s  +  2s   ,        0  <_  s  <   1 
0        ,         s  >  1 

To  obtain  the  value  of  the  interpolating  function  F(x,y)   at  the  point 

(x  ,y  ),  we  calculate  a  weighted  average  of  the  values  of  Q  (x  ,y  )   for 

*  * 
those  values  of  n  such  that   (x  ,y  )   lies  in  the  circle  C   .   In  particu- 
lar, 

F(x*,y*)  =  I     w(V(x*-xn)2  +  (y*-yn)2/rjQn(x*>y*)/l   wN(x*-xn)2  +  (7%^/^   ) 
n=l  /  n=l 

i  (7~         2     *   71 1       \  *  * 

Because   wry(x  -x)  +  (y  -y  )  /r  J     is  zero  unless   (x  ,y  )   lies  inside 

the  circle  C   ,   the  sum  can  be  taken  over  only  those  values  of  n  such 

n  J 

*  *  *  * 

that   (x  ,y  )   is  in  C   .   The  function  is  defined  for  all  points   (x  ,y  ) 

which  lie  in  the  interior  of  one  or  more  of  the  circles.   By  virtue  of  the 
weights  being  zero,  with  zero  first  partial  derivatives,  on  the  boundary  of 
the  circles,  the  interpolating  function  F  is  continuous,  with  continuous 
first  partial  derivatives. 

We  can  note  that  F  is  a  bivariate  analogue  of  the  quintic  spline  of 
deficiency  three  [1].   Unfortunately,  unlike  the  univariate  case,  the  bivar- 
iate version  above  can  fail.   It  is  not  necessarily  true  that  the  polynomials 
Q,   must  exist.   In  fact,  on  a  uniform  grid  a  boundary  point  and  its  five 
nearest  neighbors  will  lie  on  two  straight  lines,  hence  the  interpolating 
polynomial  does  not  generally  exist.   Non-existence  of  a  particular  Q,   can 
be  bypassed  by  simply  ignoring  that  circle  in  the  definition  of  F  .   This  is 


handled  computationally  by  setting   r  =  0  .   For  isolated  instances  the  effect 
may  be  minimal,  although  failure  for  very  many  points  could  severely  curtail 
the  region  of  definition  of  the  interpolating  function  and  even  exclude  some 
data  points. 

Experience  with  this  version  of  the  computational  scheme  indicates  that 
the  behavior  of  the  Q,   are  adversely  affected  by  the  presence  of  even 
moderately  steep  gradients.   This  problem  is  not  unique  to  polynomials,  but 
they  do  seem  to  be  more  sensitive  than  other  basis  functions  we  have  investi- 
gated.  In  addition,  regions  where  the  data  points  are  relatively  sparse  seem 
to  yield  poor  approximations. 

2.2  A  Generalization 

We  will  develop  a  class  of  interpolation  methods  based  on  Maude's  idea, 
which  will  contain  his  method  as  a  special  case.   We  will  begin  with  a  weight- 
ing function  associated  with  each  point,  although  in  the  application  it  will 
probably  not  be  the  starting  point. 

We  presume  that  data  points   (x,  ,y,  )   and  corresponding  function  values 

f,  ,  k  =  1 N  have  been  given.  With  each  point   (x  ,y  )   associate  a 

k  n  n 

nonnegative  weight  function  w  (x,y)  which  has  compact  support  S   ,  such 

that   (x  ,y  )   is  contained  in  the  interior  of  S   •   Let  U  (x,y)  be  a  local 
n  'n  n         n  J 

interpolating  function  which  is  defined  on  S   and  interpolates  to  the  value 

n 

f   at  each   (x,  ,y  )   in  S 

We  define  the  value  of  the  interpolating  function  F(x,y)   at  the  point 

*  * 
(x  ,y  )   to  be 

N  N 

F(x  ,y  )  -  I     wn(x  ,y  )UQ(x  ,y  )  /   J.  wn<x  »Y  )   » 

n=l  n=l 


*  * 
provided   (x  ,y  )   lies  in  the  interior  of  one  of  the  S   .   Since  pre- 

sumably   (x  ,y  )   does  not  lie  in  all  S   ,  we  may  want  to  write 


F(x  ,y  )   =   I       w  (x  ,y  )U  (x  ,y  )/  J,  w  (x  ,y  ) 
nel*  n       nel* 


*  *  * 

where  1  =  {  n  :  (x  ,y  )  e  interior  (S  )  }  . 

We  can  make  a  number  of  observations  about  the  function  F  . 

(a)  F(\»yk)  =  fk  since  if   <Vyk*  €  interior  of  S^   ,   Un<Vyk*  =  fk  » 

(b)  F  is  at  least  as  smooth  as  the  least  smooth  of  the  functions  w,   and 

*  * 
U,  ,   k  =  1,...,N  ;   (c)   F(x  ,y  )   is  locally  determined  by  the  points 

(x,  ,y,  )  e   u   S  ;  (d)  the  error  in  the  interpolation  is  no  worse  than 
k  k       .   n 
nel* 

that  of  the  poorest  interpolation  by  a  U,  ,  since  the  value  of  F  is  a 
convex  combination  of  the  U,  ,  hence  the  error  is  a  convex  combination  of 
the  local  interpolating  errors. 

It  is  trivially  observed  that  Maude's  method  is  a  special  case,  pro- 
vided the  Q,   (U,   above)  exist.   For  Maude's  method  the  error  is  of  the 

2 
form  0(h  )  ,  where  h  =  max  r,   [12].   If  one  considers  a  sequence  of  ap- 

k   k 

proximations  as  more  points  are  added,  in  a  fashion  such  that  h  ■>  0  ,  it 

is  necessary  that  the  'condition'  of  the  sets  of  nearest  neighbors  not  de- 

2. 

teriorate  too  much  in  order  to  maintain  the  0(h)   error  estimate. 


3.0  Implementation 

In  order  for  a  scheme  such  as  that  outlined  in  the  previous  section  to 
be  reasonably  easy  to  implement,  the  weight  functions,  in  connection  with 
the  shape  of  the  support  regions  S,   must  be  chosen  with  some  care.   Since 
we  wish  to  be  able  to  incorporate  smoothness,  which  requires  the  value  of 
w,   and  some  of  its  partial  derivatives  to  be  zero  on  the  boundary  of  S   , 
we  are  led  rather  naturally  to  such  uncomplicated  shapes  as  squares  and 

6 


circles.   Squares  are  associated  with  t^     and  -£..   norms  on  2-space  and 
circles  with  the  -£_  norm.   The  weight  functions  for  these  regions  can  be 
easily  built  up  from  univariate  functions,  as  follows.   Let   w(s)   denote  a 
univariate  function  with  appropriate  smoothness.   For  example 


w(s)   = 


1-|  8 
0 


8  I  <  1 

sl  >  1 


is  continuous,  while 


u(s)   = 


2    1.3 
l-3s  +  2  s 


s  |  <.  1 
s|  >  1 


has  continuous  derivatives  as  well.   Similar  functions  with  arbitrarily 
high  smoothness  can  be  generated.   Reasonable  weight  functions  for  balls 
of  radius  r,   centered  at   (x,  ,y,  )   are  given  by: 


/fr^)  +  (y-yk)  \    /  (^\)  -  (yyk)  \ 
^(x.y)-  »[ )»[ j 

wk(x,y)  =  u  L/Cx-x^)     +   (y~yk)  /rk  1 


wk(x,y)  =  w 


R)-m  ■ 


The  choice  of  region  will  affect  the  approximation  slightly.  Which  should 
be  used  is  a  matter  of  personal  preference,  in  the  author's  opinion.   The 
square  regions  have  the  distinction  that  the  weight  functions  are  piecewise 
polynomials,  which  results  in  the  function  F  being  a  rational  function  if 


the  U   are  polynomials  or  rational  functions.   While  this  is  aestheti- 
cally  more  pleasing  than  a  function  involving  square  roots  as  is  obtained 
for  circles,  the  author  believes  the  practical  difference  to  be  minimal. 
It  may  be  an  advantage  to  use  regions  which  do  not  have  flat  sides.   Some 
calculations  have  been  done  both  ways  with  neither  one  seeming  to  be  sig- 
nificantly better  than  the  other. 

We  will  generally  assume  that  the  regions  S,   are  chosen  to  contain 
a  set  number  of  nearest  neighbors,  say  m-1  ,  of  the  point   (x,  ,y,  ).   This 
will  allow  the  use  of  the  same  type  of  local  interpolating  function  for 
each  region  and  will  simplify  the  coding  of  programs  to  implement  the  pro- 
cedure.  The  value  of  m  will  greatly  influence  the  computational  effort 
required,  since  the  larger  m  is,  the  more  the  support  regions  S,   will 
overlap.   On  the  other  hand,  too  small  a  value  for  m  will  not  allow  the 
local  interpolating  function  to  mimic  f (x,y)   suitably,  and  may  cause  holes 
in  the  region  on  which  F(x,y)   is  defined,  particularly  if  the  data  is 
sparsely  located  in  some  areas.   Hence  a  balance  must  be  struck  between 
computational  efficiency  and  adequate  sampling  of  f(x,y)   along  with  cov- 
erage considerations. 

The  choice  of  local  interpolating  functions  U,   is  very  broad,  since 
essentially  the  only  restrictions  are  appropriate  smoothness  and  interpola- 
tion of  the  data.   The  obvious  choices  are  polynomials  consisting  of  m 
monomials,  chosen  in  some  fashion.   While  this  works  well  for  some  sets  of 
data,  it  has  not  yielded  the  best  approximations  overall. 

Better  results  have  been  obtained  by  using  the  optimal  approximations 
previously  mentioned  [6],  [9],  [2],  While  optimal  approximations  have  not 
received  great  attention  in  applications,  they  do  have  very  desirable  proper- 
ties and  prove  to  yield  quite  good  results,  with  one  exception,  when  applied 
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to  our  test  problems  in  the  next  section.   These  functions  are  splines  and 
thus  are  better  able  to  adapt  to  certain  data. 

When  using  optimal  local  interpolating  functions  it  is  probably  de- 
sirable to  use  weight  functions  of  no  more  smoothness  than  the  interpolating 
functions.   Thus  if  the  interpolating  function  is  the  optimal  approximation 
from  the  class   B...  ^  ,   (see  [10]  for  the  definition  of  this  space)  one  should 


1- I s  |       Is  I  <  1 


choose     w(s)   = 


are  only  continuous, 


since  functions  in  B,   , 
s  I  >   1  ' 


4.0  Numerical  Experiments 

It  was  deemed  desirable  to  test  a  number  of  methods  on  a  variety  of  sur- 
faces for  purposes  of  comparison.   Five  sets  of  data  were  fit  using  eleven 
fitting  functions,  and  one  set  of  data  was  fit  using  ten  of  these  procedures, 
the  eleventh  not  being  suitable.   We  will  describe  the  procedures  and  the  data 
sets  and  give  a  table  of  results  with  some  discussion  and  observations. 

4.1  The  Data  Sets 

Because  of  the  variety  of  surfaces  used  by  McLain  [5],  it  was  decided  to 
use  them  for  our  tests  also.   Unlike  McLain,  we  did  not  specify  uniformly  spaced 
points  on   [1,10]  x  [1,10]  .   Instead,  in  each  of  the  81  squares  [i,i+l]x[j , j+l] 
i,j  =  1,...,9   ,  a  point  was  generated  by  a  random  number  generator.   These 
points  are  shown  graphically  in  Figure  1.   In  addition,  nine  points  were  speci- 
fied at  positions  which  were  judged  a  priori  to  be  somewhat  critical  points  in 
defining  the  surface.   Thus  a  total  of  90  points  were  specified  on  surfaces 
S1-S5  on  the  square   [1,10]  x  [1,10]. 
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Figure  1:   Interpolation  Points 
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We  point  out  that  McLain  was  interested  in  contours  for  topographical 
data  and  he  describes  the  surfaces  in  topographical  terms.   We  shall  do  this 
also  as  this  leads  to  an  easy  visualization  of  the  surfaces. 

(SO   This  surface  is  a  section  of  a  sphere  and  is  described  by 


V2  2 

64  -  (x  -  5.5)  -  (y  -  5.5)    .   The  additional  points  speci- 
fied were  (cx,B)  ,  a, 3  =  1,  5.5,  10. 

(52)  This  surface  is  a  hill  rising  steeply  from  a  plain,  and  is  described 

2        2 
by  f(x,y)  =  exp(-(x-5)  -  (y-5)  )  .   The  additional  points  specified  were 

(4,5)  ,  (4.5,  4.5)  ,  (4.5,  5.5)  ,  (5,4)  ,  (5,5)  ,  (5,6)  ,  (5.5,  4.5)  ,  (5.5,  5.5), 

(6,5)  . 

(53)  This  surface  is  similar  to  S2,  a  hill  rising  less  steeply  from  a 
plain.   It  is  described  by  f(x,y)  =  exp  (-  ^-^ — *  V~5'      J.   The  additional 
points  specified  were  the  same  as  for  S2. 

(54)  This  surface  is  a  long  narrow  hill  rising  from  a  plain,  and  is  de- 

2        2 
scribed  by    f(x,y)  =  exp(-(x+y  -  11)   -  (x-y)  /10)  .   The  additional  points 

were   (3,8),  (4.5,  6.5),  (4.5,  5.5),  (5.5,  4.5),  (5.5,  5.5),  (5.5,  6.5), 

(6.5,  4.5),  (6.5,  5.5),  (8,3)  . 

(55)  This  surface  is  a  plain  and  a  plateau  separated  by  a  sharp  rise, 
almost  a  cliff.   It  is  described  by  f (x,y)  =  tanh(x  +  y  -  11)  .   The  addi- 
tional points  were  (1.10),  (2.5,  8.5),  (4,7),  (5,5),  (5.5,  5.5),  (6,6),  (7,4), 
(8.5,  2.5),  (10,1)  . 

(56)  This  surface  is  that  given  by  Ferguson  [3]  and  is  not  described 
mathematically.   For  purposes  of  better  point  spacing  the  y  coordinate  was 
multiplied  by  three.   This  was  almost  necessary  because  a  point  and  its 
nearest  five  neighbors  would  often  include  five  points  very  nearly  on  a 
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straight  line  when  the  original  data  was  used.   Multiplication  of  the  y-coordinc 
by  a  constant  is  equivalent  to  using  ellipsoidal  regions  in  methods   Fl  -  F5.   II 
has  no  effect  on  methods  F6  -  F8,  and  has  a  very  beneficial  effect  on  methods 
F9  -  F10.   The  data  is  listed  in  Table  1. 

A. 2  The  Interpolation  Functions 

The  eleven  interpolation  functions  consist  of  five  different  implemen- 
tations  of  the  method  described  in  Section  2,  three  global  methods  for 
comparison  purposes,  and  three  versions  of  McLain's  distance  weighted  least 
squares  approximation.  We  will  discuss  each  of  these  in  turn.   Several  of 
them  involve  use  of  optimal  approximations  in  certain  spaces  of  functions,  and 
references  to  more  information  on  these  topics  are  [10],  [6],  [9],  and  [2]. 

(Fl)   This  is  Maude's  method,  using  U,   as  a  second  degree  polynomial 
interpolating  to  the  point  and  its  five  nearest  neighbors.   This  was  described 
in  Section  2.1. 

(F2)   This  is  a  method  of  Section  2,  using  the  optimal  approximation 
from  the  class  of  functions  Bfl  1 ,   for  the  local  interpolating  function  for 
a  point  and  its  five  nearest  neighbors.   Six  points  per  region  was  chosen  as 
being  sufficient  for  coverage,  and  adequate  for  function  definition,  as  well 
as  comparable  computationally  to  Maude's  method.   Circular  regions  were  used 

II—  |  s  |         |  s  |  <_  1 
0  |s|  >  1 

Since  Barnhill  and  Nielson  [2]  have  shown  that  Sard  spaces  of  type  B* 
have  a  reproducing  kernel,  one  can  find  the  optimal  approximation  as  a  linear 
combination  of  basis  functions  (or  representers  of  the  point  evaluation 
functionals) ,  one  of  which  is  associated  with  each  data  point.   The  basis 
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2.0  15.0  2.0 

2.49  7.647  2.7 

2.981  0.291  2.9 

3.471  -7.062  3.0 

3.961  -14.418  3.0 

7.45  12.003  2.0 

7.35  6.012  3.0 

7.251  0.018  2.5 

7.151  -5.973  1.5 

7.051  -11.967  2.0 

10.901  9.015  1.5 

10.751  4.536  1.425 

10.602  .06  1.35 

10.453  -4.419  1.276 

10.304  -8.895  1.2 

14.055  10.509  1.0 

14.194  6.783  0.8 

14.331  3.054  1.2 

14.469  0.672  1.6 

14.607  -4.398  1.25 

15.0  12.0  0.0 

15.729  8.067  0.0 

16.457  4.134  0.20 

17.185  0.198  0.60 

17.914  -3.735  1.2 


Table     1 
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function  for  the  point      (x,  ,y   )      in     B..-    ..,      is 

\(x,y)   =    [1  +   (x-a)   -    (x-xk)+   ][1  +   (y-b)   -    (yyk>+   ]      , 


I  0   ,   t  £  0 
where    (*-)+  =  {   -r  *-s  *-he   truncated  power  function  and 

t   ,    t  >  0 


(a,b)   is  a  parameter  pair  which  in  this  case  must  satisfy  x^  >_  a  , 


n 
y  >_  b  for  n  =  1,2,...,N  .   In  general   (a,b)  need  not  satisfy  this 


condition,  but  we  have  simplified  K,   by  requiring  the  condition. 

We  can  note  that  optimal  approximations  from  B.-  -,   are  piecewise 
bilinear  functions  with  knots  along  lines  through  the  data  points  and  paral- 
lel to  the  axes. 

(F3)   This  is  a  method  of  Section  2,  similar  to  function  F2,  except  that 


U,   is  taken  to  be  the  optimal  approximation  from  B,^  ?,  ,  and  the  weight 
function  used    w(s)   = 


l-3s2  -4-  2  Is  |        Is  I  <  1 


0  s  >  1 


Again,  the  nearest  five  neighbors  in  a  circular  region  were  used.   The  basis 
function  associated  with  the  point  (xu»y\c)      *n  this  case  is 

^(x.y)  =  g(a,xk,x)g(b,yk,y) 

(x^a)  (x-a)      (x-a)    (x~\)+ 
where     g(a,xfc,x)  -  1  +  (x^  -  a)  (x-a)  +  jj "  — Jj —  +  — J] i 

where  the  notation  and  parameters  are  as  described  before. 

(F4)   This  approximation  is  the  same  as  F3  except  that  the  nine  nearest 
neighbors  in  a  circular  region  were  used  to  define  the  local  interpolating 
function.   This  is  to  test  the  effect  of  including  more  points  in  the  local 
approximation. 
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(F5)   This  is  a  method  of  Section  2,  where  the  U,   is  taken  to  be  the 

optimal  approximation  from  the  space  of  functions  T  '   described  by 
Mansfield  [6].   Again,  a  circular  region  containing  a  point  and  its  five 

1-  |  s  |   ,   |  s  |  <_  1 
nearest  neighbors  was  used,  with   w(s)  =  {  .   This 

0    ,   |s|  >  1 

space  is  constrained  so  as  to  give  an  exact  approximation  if  the  approxi- 
mated function  is  linear.   Our  version  used  the  basis  functions  given  by 
Nielson  [9]   and  are  of  the  form 

V     <  N     (X~V+   +   ^k^   + 

Vx,y)  "  —3! +  3I-  + 

[(xk  "  x>+-  <a  "  x>+  "  (\  "  a)+  ]  [  <  V>+  -  (b-y)+  -  (Vb)J  . 

where  a  and  b  are  parameters.   The  optimal  approximation  is  a  linear 
combination  of  l,x,y,  and  the  K,  . 

(F6)   This  is  the  global  optimal  approximation  from  B...  - ,  . 

(F7)   This  is  the  global  optimal  approximation  from  B.-  «,  . 

(F8)   This  is  the  global  optimal  approximation  from  T1'1  . 

(F9)   This  is  McLain's  method  M10  which  he  describes  as  being  accurate 

with  a  satisfactory  amount  of  computation.   The  weight  attached  to  a  data 

2   2  2        2        2 

point   (^yj)   is  w(d)  =  exp(d  )/d    where  d  =  (x-x,  )  +  (y_yfc)   is  the 

distance  squared  from  the  point   (x,y)   at  which  the  function  value  is  desired, 

The  fitting  function  is  a  second  degree  polynomial. 

(F10)   This  is  the  same  as  F9  except  that  distance  is  changed.   Here 

2         2        2 
d  =  [(x-x.)  +  (y-y,)  ]/10  .   This  is  equivalent  to  shrinking  the  coordinates 

by  a  factor  YlO  ,  and  is  to  test  the  effect  of  a  different  coordinate  system. 
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d2  =  10 


2        2 


(Fll)   This  is  the  same  as  F9  except  that  distance  is  changed.   Here 

This  is  equivalent  to  expanding  the  coordi- 
nates by  a  factor  /10  and  is  to  test  the  effect  of  a  different  coordinate 
system. 

4. 3  Comparisons 

For  surfaces  S1-S5  and  all  eleven  interpolating  functions  the  devi- 
ation at  the  361  points   (x,y)  ,  x,y  =  1(.5)10  was  generated.   The  maximum 
deviation,  the  mean  deviation,  and  the  root-mean-square  deviations  are  listed 
in  Table  2.   The  number  in  parenthesis  below  each  function  is  the  approximate 
calculation  time  in  seconds  for  all  calculated  data  points.   The  calculations 
were  done  on  the  IBM  360/67  at  the  Naval  Postgraduate  School.   These  times 
are  given  only  as  an  indication  of  the  number  of  computations  required  for 
the  various  methods. 

Generally  we  can  observe  several  things  from  the  table: 

(i)    Maude's  method  (Fl)  does  not  compare  favorably  with  any  of  the 
other  methods,  except  when  applied  to  surface  SI  . 

(ii)   Surface  SI  seems  to  be  relatively  difficult  to  fit,  although 
Maude's  method   (Fl)   and  McLain's  method   (F9)  work  very  well. 

(iii)   The  local  weighted  optimal  approximations   (F2-F5)   often  are  more 
accurate  than  the  corresponding  global  versions   (F6-F8)  .   A  notable  excep- 
tion is  the  application  to  surface  Si  . 

(iv)   McLain's  method   (F9)   does  very  well  on  all  these  surfaces  and 
is  often  the  most  accurate.   The  suspicion  that  the  distance  weighting  func- 
tion is  rather  finely  tuned  for  this  particular  coordinate  system  is  rein- 
forced by  poorer  performance  of  the  expanded  and  contracted  length  versions 
of  it   (F10-F11) .   The  amount  of  computation  required  is  considerably  greater 
for  the  method  than  for  any  other  methods  considered  here,  even  the  global 
methods. 
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Surface 

SI 

S2 

S3 

S4 

S5 

Fitting  Function 

Fl  :  Maude 
(34) 

.06562 
.00390 
.00791 

1.390 
.06882 
.2106 

1.198 
.03964 
.1262 

3.798 
.09371 
.3789 

.7438 

.04440 

.10468 

F2  :  Bj.    .j  local 
(51) 

1.768 
.1378 
.2456 

.1704 

.00943 

.02661 

.1520 

.01180 

.02429 

.3393 

.01732 

.04202 

.4926 

.04828 

.09124 

F3  :  Br2)2l  local 
(63) 

1.284 
.1049 
.1859 

.1390 

.00726 

.01993 

.08185 
.00832 
.01483 

.2701 

.01446 

.03676 

1.111 
.03876 
.09576 

F4  :  Br2>21  local 
(143) 

.3472 

.03970 

.06454 

.1315 

.00727 

.01953 

.05723 
.00593 
.01073 

.2585 

.01758 

.03831 

1.178 
.04932 
.11336 

F5  :  T1,1  local 
(60) 

.1068 

.01667 

.02477 

.1740 

.00884 

.02149 

.1143 

.01035 

.01797 

.2398 

.01493 

.03280 

.4656 

.03872 

.07300 

F6  :  Bn  _,  global 
(128) 

.2209 

.02724 

.05102 

.2829 

.05706 

.08218 

.1514 

.02234 

.03416 

.3547 

.05453 

.08126 

.5098 

.10322 

.14568 

F7  :  Br2>2l  global 
(154) 

.09152 
.00832 
.01465 

.4838 

.07447 

.1086 

.09318 
.01012 
.01719 

1.335 
.2005 
.3127 

.8164 

.10104 

.15686 

F8  :  T1'1  global 
(138) 

.1222 

.01355 

.02394 

.2693 

.05135 

.07601 

.1457 

.02001 

.03150 

.3549 

.05106 

.07797 

.4868 

.09632 

.13595 

F9  :  McLain 
(372) 

.05223 
.00354 
.00827 

.1372 

.00824 

.02090 

.04894 
.00530 
.00873 

.1803 

.01124 

.02546 

.4164 

.03006 

.05960 

F10:  McLain,  d2  +   d  /10 
(465) 

.04343 
.00621 
.00982 

.2158 

.02196 

.04363 

.08728 
.02174 
.03117 

.3894 

.02549 

.05174 

.2738 

.05374 

.07886 

2      2 
Fll:  McLain,  d  *•   ^Qd 

(87) 

5.735 
.03785 
.3169 

1.304 
.01222 
.07502 

.1254 

.00982 

.02010 

.4145 

.01539 

.03949 

4.118 
.04584 
.2318 

Table  2:   Maximum,  Mean  and  RMS  deviations 
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(v)   While  functions  in  B.   .,   and  T  *    are  of  similar  smoothness 
1,1       l-L.J-1 
(only  continuous),   T  !    gives  considerably  more  flexibility  which  shows 

up  when  fitting  surface  SI.   Results  are  comparable  for  surfaces   S2-S5. 


(vi)  Changing  the  number  of  neighbors  included  in  a  region  seems  to 
generally  have  only  a  small  effect   (F3  vs.  F4) .   The  one  exception,  again, 
is  on  surface  Si  . 

For  surface  S6  we  could  not  compare  exact  errors  since  only  data 
points  were  given.   The  data  is  apparently  difficult  to  fit  since  it  in- 
volves three  relative  maxima  with  a  saddle  point  between  two  of  them  and 
a  shallow  valley  between  two  of  them,  the  entire  surface  only  having  25 
points  given  on  it.   In  this  situation  it  is  not  surprising  that  interpo- 
lating surfaces  may  involve  some  overshoot.   We  list  the  minimum  and  maxi- 
mum observed  values  for  functions  F1-F10  in  Table  3.   Again,  the  flexi- 
bility of  the  space  T  '   is  noted,  and  the  fine  tuning  of  the  distance 
function  for  F9  ,  yielding  a  good  result  in  F10. 
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Function 
Fl 
F2 
F3 
F4 
F5 
F6 
F7 
F8 
F9 
F10 


Extreme  Values 

-6.43, 

7.18 

0   , 

4.65 

0   , 

5.52 

0   , 

6.02 

0  , 

3.00 

0   , 

3.77 

0   , 

5.05 

0  , 

3.05 

-26.31, 

25.83 

0   , 

3.30 

Table  3 
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5.0   Conclusion 

The  main  conclusion  which  can  be  made  is  that  the  method  of  Section  2 
yields  good  approximations  provided  the  local  interpolating  functions  are 
reasonably  accurate.   Six  points  per  region  appears  to  give  sufficient  func- 
tion definition  and  coverage  at  a  reasonable  computational  load.   Any  of  the 
optimal  approximations  seem  to  be  satisfactory  as  local  interpolating  func- 
tions, and  the  use  of  them  is  strongly  recommended  over  polynomials. 

The  use  of  McLain's  method  (F9)  or  a  variation  of  it  gives  very  good 
results  if  the  computational  load  is  not  burdensome.   Some  tuning  of  the 
measure  of  distance  may  be  necessary  and  is  a  definite  disadvantage  of  the 
method. 

It  appears  that  on  certain  surfaces,  some  methods  will  do  poorly  while 
on  others  the  situation  could  be  reversed.   In  terms  of  which  of  the  methods 
discussed  can  ultimately  be  tuned  to  a  given  set  of  data,  the  author  feels 
McLain's  method  has  a  good  capability.   Unfortunately  in  case  only  isolated 
points  are  known  the  procedure  for  this  tuning  is  not  known,  and  it  is  proba- 
bly safer  to  use  a  method  which  adapts  naturally  to  changes  in  coordinate 
systems  and  possible  varying  sparseness  of  data  points,  such  as  the  methods 
of  Section  2,  where  changing  the  obvious  parameter  (the  number  of  neighbors 
included  in  a  region)  seems  to  have  a  fairly  small  influence. 
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